Brian E. J. Rose, University at Albany
You really should be looking at The Climate Laboratory book by Brian Rose, where all the same content (and more!) is kept up to date.
Here you are likely to find broken links and broken code.
This document uses the interactive Jupyter notebook format. The notes can be accessed in several different ways:
github at https://github.com/brian-rose/ClimateModeling_coursewareAlso here is a legacy version from 2015.
Many of these notes make use of the climlab package, available at https://github.com/brian-rose/climlab
In [1]:
# Ensure compatibility with Python 2 and 3
from __future__ import print_function, division
Recall from last lecture that the bulk formula for surface evaporation (latent heat flux) is
$$ \text{LE} = L ~\rho ~ C_D ~ U \left( q_s - q_a \right) $$which we approximated in terms of temperatures for a wet surface as
$$ \text{LE} \approx L ~\rho ~ C_D ~ U \left( (1-r) ~ q_s^* + r \frac{\partial q^*}{\partial T} \left( T_s - T_a \right) \right) $$The drag coefficient $C_D$ determines the flux for a given set of temperatures, relative humidity, and wind speed.
Now suppose that the drag coefficient is reduced by a factor of two (for evaporation only, not for sensible heat flux). i.e. all else being equal, there will be half as much evaporation.
Reasoning through the effects of this perturbation (and calculating the effects in models) will give us some insight into several different roles played by water in the climate system.
What is the effect of the reduced evaporation efficiency on surface temperature?
We can use climlab to construct a model for the zonal-average climate. The model will be on a pressure-latitude grid. It will include the following processes:
This model basically draws together all the process models we have developed throughout the course, and adds the surface flux parameterizations.
Note that since we are using explicit surface flux parameterizations, we will now use the convective adjustment only on the atmospheric air temperatures. Previous our adjustment has also modified the surface temperature, which was implicitly taking account of the turbulent heat fluxes.
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import climlab
from climlab import constants as const
In [3]:
def inferred_heat_transport( energy_in, lat_deg ):
'''Returns the inferred heat transport (in PW) by integrating the net energy imbalance from pole to pole.'''
from scipy import integrate
from climlab import constants as const
lat_rad = np.deg2rad( lat_deg )
return ( 1E-15 * 2 * np.math.pi * const.a**2 *
integrate.cumtrapz( np.cos(lat_rad)*energy_in,
x=lat_rad, initial=0. ) )
In [4]:
# A two-dimensional domain
num_lev = 50
state = climlab.column_state(num_lev=num_lev, num_lat=60, water_depth=10.)
lev = state.Tatm.domain.axes['lev'].points
Here we specify cloud properties. The combination of the two cloud layers defined below were found to reproduce the global, annual mean energy balance in a single-column model.
We will specify the same clouds everywhere for simplicity. A more thorough investigation would incorporate some meridional variations in cloud properties.
In [5]:
# Define two types of cloud, high and low
cldfrac = np.zeros_like(state.Tatm)
r_liq = np.zeros_like(state.Tatm)
r_ice = np.zeros_like(state.Tatm)
clwp = np.zeros_like(state.Tatm)
ciwp = np.zeros_like(state.Tatm)
# indices
high = 10 # corresponds to 210 hPa
low = 40 # corresponds to 810 hPa
# A high, thin ice layer (cirrus cloud)
r_ice[:,high] = 14. # Cloud ice crystal effective radius (microns)
ciwp[:,high] = 10. # in-cloud ice water path (g/m2)
cldfrac[:,high] = 0.322
# A low, thick, water cloud layer (stratus)
r_liq[:,low] = 14. # Cloud water drop effective radius (microns)
clwp[:,low] = 100. # in-cloud liquid water path (g/m2)
cldfrac[:,low] = 0.21
# wrap everything up in a dictionary
mycloud = {'cldfrac': cldfrac,
'ciwp': ciwp,
'clwp': clwp,
'r_ice': r_ice,
'r_liq': r_liq}
In [6]:
plt.plot(cldfrac[0,:], lev)
plt.gca().invert_yaxis()
plt.ylabel('Pressure (hPa)')
plt.xlabel('Cloud fraction')
plt.title('Prescribed cloud fraction in the column model')
plt.show()
In [7]:
# The top-level model
model = climlab.TimeDependentProcess(state=state, name='Radiative-Convective-Diffusive Model')
# Specified relative humidity distribution
h2o = climlab.radiation.ManabeWaterVapor(state=state)
# Hard convective adjustment for ATMOSPHERE ONLY (not surface)
conv = climlab.convection.ConvectiveAdjustment(state={'Tatm':model.state['Tatm']},
adj_lapse_rate=6.5,
**model.param)
# Annual mean insolation as a function of latitude and time of year
sun = climlab.radiation.DailyInsolation(domains=model.Ts.domain)
# Couple the radiation to insolation and water vapor processes
rad = climlab.radiation.RRTMG(state=state,
specific_humidity=h2o.q,
albedo=0.125,
insolation=sun.insolation,
coszen=sun.coszen,
**mycloud)
model.add_subprocess('Radiation', rad)
model.add_subprocess('Insolation', sun)
model.add_subprocess('WaterVapor', h2o)
model.add_subprocess('Convection', conv)
print( model)
Here we add a diffusive heat transport process. The climlab code is set up to handle meridional diffusion level-by-level with a constant coefficient.
In [8]:
from climlab.dynamics import MeridionalDiffusion
# thermal diffusivity in W/m**2/degC
D = 0.04
# meridional diffusivity in m**2/s
K = D / model.Tatm.domain.heat_capacity[0] * const.a**2
d = MeridionalDiffusion(state={'Tatm': model.state['Tatm']},
K=K, **model.param)
model.add_subprocess('Diffusion', d)
Now we will add the surface heat flux processes. We have not used these before.
Note that the drag coefficient $C_D$ is passed as an input argument when we create the process. It is also stored as an attribute of the process and can be modified (see below).
The bulk formulas depend on a wind speed $U$. In this model, $U$ is specified as a constant. In a model with more complete dynamics, $U$ would be interactively calculated from the equations of motion.
In [9]:
# Add surface heat fluxes
shf = climlab.surface.SensibleHeatFlux(state=model.state, Cd=0.5E-3)
lhf = climlab.surface.LatentHeatFlux(state=model.state, Cd=0.5E-3)
# set the water vapor input field for LHF
lhf.q = h2o.q
model.add_subprocess('SHF', shf)
model.add_subprocess('LHF', lhf)
In [10]:
print( model)
In [11]:
model.integrate_years(4.)
In [12]:
# One more year to get annual-mean diagnostics
model.integrate_years(1.)
In [13]:
ticks = [-90, -60, -30, 0, 30, 60, 90]
fig, axes = plt.subplots(2,2,figsize=(14,10))
ax = axes[0,0]
ax.plot(model.lat, model.timeave['Ts'])
ax.set_title('Surface temperature (reference)')
ax.set_ylabel('K')
ax2 = axes[0,1]
field = (model.timeave['Tatm']).transpose()
cax = ax2.contourf(model.lat, model.lev, field)
ax2.invert_yaxis()
fig.colorbar(cax, ax=ax2)
ax2.set_title('Atmospheric temperature (reference)');
ax2.set_ylabel('hPa')
ax3 = axes[1,0]
ax3.plot(model.lat, model.timeave['LHF'], label='LHF')
ax3.plot(model.lat, model.timeave['SHF'], label='SHF')
ax3.set_title('Surface heat flux (reference)')
ax3.set_ylabel('W/m2')
ax3.legend();
ax4 = axes[1,1]
Rtoa = np.squeeze(model.timeave['ASR'] - model.timeave['OLR'])
ax4.plot(model.lat, inferred_heat_transport(Rtoa, model.lat))
ax4.set_title('Meridional heat transport (reference)');
ax4.set_ylabel('PW')
for ax in axes.flatten():
ax.set_xlim(-90,90); ax.set_xticks(ticks)
ax.set_xlabel('Latitude'); ax.grid();
In [14]:
model2 = climlab.process_like(model)
model2.subprocess['LHF'].Cd *= 0.5
model2.integrate_years(4.)
model2.integrate_years(1.)
In [15]:
fig, axes = plt.subplots(2,2,figsize=(14,10))
ax = axes[0,0]
ax.plot(model.lat, model2.timeave['Ts'] - model.timeave['Ts'])
ax.set_title('Surface temperature anomaly')
ax.set_ylabel('K')
ax2 = axes[0,1]
field = (model2.timeave['Tatm'] - model.timeave['Tatm']).transpose()
cax = ax2.contourf(model.lat, model.lev, field)
ax2.invert_yaxis()
fig.colorbar(cax, ax=ax2)
ax2.set_title('Atmospheric temperature anomaly');
ax2.set_ylabel('hPa')
ax3 = axes[1,0]
for field in ['LHF','SHF']:
ax3.plot(model2.lat, model2.timeave[field] - model.timeave[field], label=field)
ax3.set_title('Surface heat flux anomalies')
ax3.set_ylabel('W/m2')
ax3.legend();
ax4 = axes[1,1]
Rtoa = np.squeeze(model.timeave['ASR'] - model.timeave['OLR'])
Rtoa2 = np.squeeze(model2.timeave['ASR'] - model2.timeave['OLR'])
ax4.plot(model.lat, inferred_heat_transport(Rtoa2-Rtoa, model.lat))
ax4.set_title('Meridional heat transport anomaly');
ax4.set_ylabel('PW')
for ax in axes.flatten():
ax.set_xlim(-90,90); ax.set_xticks(ticks)
ax.set_xlabel('Latitude'); ax.grid();
print ('The global mean surface temperature anomaly is %0.2f K.'
%np.average(model2.timeave['Ts'] - model.timeave['Ts'],
weights=np.cos(np.deg2rad(model.lat)), axis=0) )
This model predicts the following:
Basically, this model predicts that by inhibiting evaporation in the tropics, we force the tropical surface to warm and the tropical atmosphere to cool. This cooling signal is then communicated globally by atmospheric heat transport. The result is small positive global surface temperature anomaly.
We could list many things, but as we will see below, two key climate components that are not included in this model are
We will compare this result to an analogous experiment in a GCM.
The model is the familiar CESM but in simplified "aquaplanet" setup. The surface is completely covered by a shallow slab ocean.
This model setup (with CAM4 model physics) is described in detail in this paper:
Here we will compare a control simulation with a perturbation simulation in which we have once again reduced the drag coefficient by a factor of 2.
In [16]:
# Load the climatologies from the CAM4 aquaplanet runs
datapath = "http://ramadda.atmos.albany.edu:8080/repository/opendap/latest/Top/Users/BrianRose/CESM_runs/"
endstr = "/entry.das"
ctrl = xr.open_dataset(datapath + 'aquaplanet_som/QAqu_ctrl.cam.h0.clim.nc' + endstr, decode_times=False).mean(dim='time')
halfEvap = xr.open_dataset(datapath + 'aquaplanet_som/QAqu_halfEvap.cam.h0.clim.nc' + endstr, decode_times=False).mean(dim='time')
In [17]:
lat = ctrl.lat
lon = ctrl.lon
lev = ctrl.lev
In [18]:
TS_anom = halfEvap.TS - ctrl.TS
Tatm_anom = halfEvap['T'] - ctrl['T']
In [19]:
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(14,5))
ax1.plot(lat, TS_anom.mean(dim='lon')); ax1.set_title('Surface temperature anomaly')
cax2 = ax2.contourf(lat, lev, Tatm_anom.mean(dim='lon'), levels=np.arange(-7, 8., 2.), cmap='seismic')
ax2.invert_yaxis(); fig.colorbar(cax2,ax=ax2); ax2.set_title('Atmospheric temperature anomaly');
for ax in (ax1, ax2):
ax.set_xlim(-90,90); ax.set_xticks(ticks); ax.grid();
print ('The global mean surface temperature anomaly is %0.2f K.' %((TS_anom*ctrl.gw).mean(dim=('lat','lon'))/ctrl.gw.mean(dim='lat')))
In this model, reducing the evaporation efficiency leads to a much warmer climate. The largest warming occurs in mid-latitudes. The warming is not limited to the surface but in fact extends deeply through the troposphere.
Both the spatial pattern and the magnitude of the warming are completely different than what our much simpler model predicted.
Why?
In [20]:
energy_budget = {}
for name, run in zip(['ctrl','halfEvap'],[ctrl,halfEvap]):
budget = xr.Dataset()
# TOA radiation
budget['OLR'] = run.FLNT
budget['OLR_clr'] = run.FLNTC
budget['ASR'] = run.FSNT
budget['ASR_clr'] = run.FSNTC
budget['Rtoa'] = budget.ASR - budget.OLR # net downwelling radiation
# surface fluxes (all positive UP)
budget['LHF'] = run.LHFLX
budget['SHF'] = run.SHFLX
budget['LWsfc'] = run.FLNS
budget['LWsfc_clr'] = run.FLNSC
budget['SWsfc'] = -run.FSNS
budget['SWsfc_clr'] = -run.FSNSC
budget['SnowFlux'] = ((run.PRECSC+run.PRECSL)
*const.rho_w*const.Lhfus)
# net upward radiation from surface
budget['SfcNetRad'] = budget['LWsfc'] + budget['SWsfc']
budget['SfcNetRad_clr'] = budget['LWsfc_clr'] + budget['SWsfc_clr']
# net upward surface heat flux
budget['SfcNet'] = (budget['SfcNetRad'] + budget['LHF'] +
budget['SHF'] + budget['SnowFlux'])
# net heat flux in to atmosphere
budget['Fatmin'] = budget['Rtoa'] + budget['SfcNet']
# hydrological cycle
budget['Evap'] = run['QFLX'] # kg/m2/s or mm/s
budget['Precip'] = (run['PRECC']+run['PRECL'][:])*const.rho_w # kg/m2/s or mm/s
budget['EminusP'] = budget.Evap - budget.Precip # kg/m2/s or mm/s
energy_budget[name] = budget
In [21]:
# Here we take advantage of xarray!
# We can simply subtract the two xarray.Dataset objects
# to get anomalies for every term
# And also take the zonal averages for all anomaly fields in one line of code
anom = energy_budget['halfEvap'] - energy_budget['ctrl']
zonanom = anom.mean(dim='lon')
In [22]:
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(14,5))
ax1.plot(lat, zonanom.ASR, color='b', label='ASR')
ax1.plot(lat, zonanom.OLR, color='r', label='OLR')
ax1.plot(lat, zonanom.ASR_clr, color='b', linestyle='--')
ax1.plot(lat, zonanom.OLR_clr, color='r', linestyle='--')
ax1.set_title('TOA radiation anomalies')
ax2.plot(lat, zonanom.SWsfc, color='b', label='SW')
ax2.plot(lat, zonanom.SWsfc_clr, color='b', linestyle='--')
ax2.plot(lat, zonanom.LWsfc, color='g', label='LW')
ax2.plot(lat, zonanom.LWsfc_clr, color='g', linestyle='--')
ax2.plot(lat, zonanom.LHF, color='r', label='LHF')
ax2.plot(lat, zonanom.SHF, color='c', label='SHF')
ax2.plot(lat, zonanom.SfcNet, color='m', label='Net')
ax2.set_title('Surface energy budget anomalies')
for ax in [ax1, ax2]:
ax.set_ylabel('W/m2'); ax.set_xlabel('Latitude')
ax.set_xlim(-90,90); ax.set_xticks(ticks);
ax.legend(); ax.grid();
Dashed lines are clear-sky radiation anomalies.
Looking at the TOA budget:
This is very suggestive of an important role for low-level cloud changes. [Why?]
From the surface budget:
In [23]:
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(12,5))
RH = (halfEvap.RELHUM - ctrl.RELHUM).mean(dim='lon'); CLOUD = (halfEvap.CLOUD - ctrl.CLOUD).mean(dim='lon')
contours = np.arange(-15, 16., 2.)
cax1 = ax1.contourf(lat, lev, RH, levels=contours, cmap='seismic'); fig.colorbar(cax1, ax=ax1); ax1.set_title('Relative Humidity (%)')
cax2 = ax2.contourf(lat, lev, 100*CLOUD, levels=contours, cmap='seismic'); ax2.set_title('Cloud fraction (%)')
for ax in [ax1, ax2]:
ax.invert_yaxis(); ax.set_xlim(-90,90); ax.set_xticks(ticks);
In [24]:
HT = {}
HT['total'] = inferred_heat_transport(anom.Rtoa.mean(dim='lon'), lat)
HT['atm'] = inferred_heat_transport(anom.Fatmin.mean(dim='lon'), lat)
HT['latent'] = inferred_heat_transport(anom.EminusP.mean(dim='lon') * const.Lhvap, lat)
HT['dse'] = HT['atm'] - HT['latent']
In [25]:
fig, ax = plt.subplots()
ax.plot(lat, HT['total'], 'k-', label='total', linewidth=2)
ax.plot(lat, HT['dse'], 'b', label='dry')
ax.plot(lat, HT['latent'], 'r', label='latent')
ax.set_xlim(-90,90); ax.set_xticks(ticks); ax.grid()
ax.legend(loc='upper left'); ax.set_ylabel('PW'); ax.set_xlabel('Latitude')
Out[25]:
We have forced a climate change NOT by adding any kind of radiative forcing, but just by changing the efficiency of evaporation at the sea surface.
The climate system then find a new equilibrium in which the radiative fluxes, surface temperature, air-sea temperature difference, boundary layer relative humidity, and wind speeds all change simultaneously.
Reasoning our way through such a problem from first principles in practically impossible. This is particularly true because in this example, the dominant driver of the climate change is an increase in SW absorption due to a substantial decrease in low-level clouds across the subtropics and mid-latitudes.
A comprehensive theory to explain these cloud changes does not yet exist. Understanding changes in low-level cloudiness under climate change is enormously important -- because these clouds, which have an unambiguous cooling effect, are a key determinant of climate sensitivity. There is lots of work left to do.
Water is intimately involved in just about every aspect of the planetary energy budget. Here we have highlighted the role of water in:
In [26]:
%load_ext version_information
%version_information numpy, scipy, matplotlib, xarray, climlab
Out[26]:
The author of this notebook is Brian E. J. Rose, University at Albany.
It was developed in support of ATM 623: Climate Modeling, a graduate-level course in the Department of Atmospheric and Envionmental Sciences
Development of these notes and the climlab software is partially supported by the National Science Foundation under award AGS-1455071 to Brian Rose. Any opinions, findings, conclusions or recommendations expressed here are mine and do not necessarily reflect the views of the National Science Foundation.
In [ ]: